//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Uniform_Equal_3D.cc
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Uniform_Equal_3D member definitions.
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <cmath>
#include <cstdio>

#include "angle/Uniform_Equal_3D.hh"
#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"

namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
/*!
 * \brief Constructor.
 */
Uniform_Equal_3D::Uniform_Equal_3D( int    sn_order,
		                            double norm,
		                            bool   uniform_xi )
: Base(sn_order, norm)
{
	using util::soft_equiv;
    using constants::pi;
    using constants::four_pi;


	Require (d_sn_order >= 2);
	Require (d_norm > 0.0);
	Require (d_sn_order % 2 == 0);

	// size angles 2D/3D
	d_angles.resize(3);

	// set pointers to mu, eta, xi
	Base::Vec_Dbl *mu  = &d_angles[MU];
	Base::Vec_Dbl *eta = &d_angles[ETA];
	Base::Vec_Dbl *xi =  &d_angles[XI];

	// number of polar angles
	int num_xi = d_sn_order;

	// number of azimuthal angles
	int num_omega = 2 * num_xi;

	// this should be satisfied (even d_sn_order and num_omega = num_xi*2)
	Require (num_omega % 4 == 0);

	// number of octants
	int num_oct = 8;

	// total number of angles in the quadrature
	int num_ang = num_xi * num_omega;

	// number of angles we will be using (2D/3D-dependent)
	int num_ang_store = num_ang / 8 * num_oct;

	// common values below
	int num_omega_4 = num_omega/4;
	int num_xi_2    = num_xi/2;
	double two_num_xi = 2.0 / num_xi;

	// temporary arrays to hold the positive octant values
	double tmp_xi[num_xi_2+1]; // +1 to handle order=2 case
	double tmp_mu[num_xi_2][num_omega_4];
	double tmp_eta[num_xi_2][num_omega_4];

	// uniform xi density
	if ( uniform_xi )
	{
		// uniform values of the cosine of the polar angle
		for ( int i = 0; i < num_xi_2+1; i++ ) // +1 to handle order=2 case
			tmp_xi[i] = two_num_xi*(i+0.5+0.5*num_xi)-1.0;
	}
	else
	{
		// uniform values of the polar angle
		// NOT IMPLEMENTED
		throw util::assertion(
				"Option for uniform polar angle density not implemented.");
	}

	// width of polar angle cosine mesh
	double del_xi = tmp_xi[1]-tmp_xi[0]; // size of tmp_xi >= 2 from above +1's

	// temporary azimuth cosine and sine arrays
	tmp_mu[num_omega_4][num_xi_2];
	tmp_eta[num_omega_4][num_xi_2];

	double coef = 4.0*pi/(num_ang*del_xi);
	double tmp  = 0;
	double cos_omega = 0;
	double sin_omega = 0;

	// compute mu and eta
	for ( int j = 0; j < num_omega_4; j++ )
	{
		tmp       = coef*(j+0.5);
	    cos_omega = cos(tmp);
	    sin_omega = sin(tmp);
	    for ( int i = 0; i < num_xi_2; i++ )
	    {
	        tmp_mu[j][i]  = cos_omega * sqrt(1.0-tmp_xi[i]*tmp_xi[i]);
	        tmp_eta[j][i] = sin_omega * sqrt(1.0-tmp_xi[i]*tmp_xi[i]);
	    }
	}

	// Force the direction vectors to be the correct length.
	mu->resize(num_ang_store);
	eta->resize(num_ang_store);
	xi->resize(num_ang_store);
	d_wt.resize(num_ang_store);

	// All angles have a constant weight
	double wt = four_pi / num_ang_store;

	// Octant directions
	double s[8][3] = { {1,1,1},{-1,1,1},{1,-1,1},{-1,-1,1},
			           {1,1,-1},{-1,1,-1},{1,-1,-1},{-1,-1,-1} };

	// Fill the quadrature angle and weight arrays
	int k = 0;
	for ( int o = 0; o < num_oct; o++ )
	{
	    for ( int i = 0; i < num_xi_2; i++ )
	    {
	        for ( int j = 0; j < num_omega_4; j++ )
	        {
	            (*mu)[k]  	= s[o][MU]*tmp_mu[j][i];
	            (*eta)[k] 	= s[o][ETA]*tmp_eta[j][i];
	            (*xi)[k]  	= s[o][XI]*tmp_xi[i];
	            d_wt[k] 	= wt;
	            ++k;
	        }
	    }
	}
	Ensure ( k == num_ang_store );

	// Normalize the quadrature set
	double wsum = 0.0;
	for (size_t n = 0; n <= num_ang_store - 1; ++n)
		wsum = wsum + d_wt[n];

	for (size_t n = 0; n <= num_ang_store - 1; ++n)
		d_wt[n] = d_wt[n] * (norm / wsum);

#ifdef DEBUG
	display();
#endif

	Ensure (d_angles[MU].size() == num_ang_store);
	Ensure (d_angles[ETA].size() == num_ang_store);
	Ensure (d_angles[XI].size() == num_ang_store);
}

/*!
 * \brief Constructor.
 */
Uniform_Equal_2D::Uniform_Equal_2D( int    sn_order,
		                            double norm,
		                            bool   uniform_xi )
: Base(sn_order, norm)
{
	using util::soft_equiv;
    using constants::pi;
    using constants::four_pi;

	Require (d_sn_order >= 2);
	Require (d_norm > 0.0);
	Require (d_sn_order % 2 == 0);

	// size angles 2D/3D
	d_angles.resize(2);

	// set pointers to mu, eta, xi
	Base::Vec_Dbl *mu  = &d_angles[MU];
	Base::Vec_Dbl *eta = &d_angles[ETA];

	// number of polar angles
	int num_xi = d_sn_order;

	// number of azimuthal angles
	int num_omega = 2 * num_xi;

	// this should be satisfied (even d_sn_order and num_omega = num_xi*2)
	Require (num_omega % 4 == 0);

	// number of octants
	int num_oct = 4;

	// total number of angles in the quadrature
	int num_ang = num_xi * num_omega;

	// number of angles we will be using (2D/3D-dependent)
	int num_ang_store = num_ang / 8 * num_oct;

	// common values below
	int num_omega_4 = num_omega/4;
	int num_xi_2    = num_xi/2;
	double two_num_xi = 2.0 / num_xi;

	// temporary arrays to hold the positive octant values
	double tmp_xi[num_xi_2+1]; // +1 to handle order=2 case
	double tmp_mu[num_xi_2][num_omega_4];
	double tmp_eta[num_xi_2][num_omega_4];

	// uniform xi density
	if ( uniform_xi )
	{
		// uniform values of the cosine of the polar angle
		for ( int i = 0; i < num_xi_2+1; i++ ) // +1 to handle order=2 case
			tmp_xi[i] = two_num_xi*(i+0.5+0.5*num_xi)-1.0;
	}
	else
	{
		// uniform values of the polar angle
		// NOT IMPLEMENTED
		throw util::assertion(
				"Option for uniform polar angle density not implemented.");
	}

	// width of polar angle cosine mesh
	double del_xi = tmp_xi[1]-tmp_xi[0]; // size of tmp_xi >= 2 from above +1's

	// temporary azimuth cosine and sine arrays
	tmp_mu[num_omega_4][num_xi_2];
	tmp_eta[num_omega_4][num_xi_2];

	double coef = 4.0*pi/(num_ang*del_xi);
	double tmp  = 0;
	double cos_omega = 0;
	double sin_omega = 0;

	// compute mu and eta
	for ( int j = 0; j < num_omega_4; j++ )
	{
		tmp       = coef*(j+0.5);
	    cos_omega = cos(tmp);
	    sin_omega = sin(tmp);
	    for ( int i = 0; i < num_xi_2; i++ )
	    {
	        tmp_mu[j][i]  = cos_omega * sqrt(1.0-tmp_xi[i]*tmp_xi[i]);
	        tmp_eta[j][i] = sin_omega * sqrt(1.0-tmp_xi[i]*tmp_xi[i]);
	    }
	}

	// Force the direction vectors to be the correct length.
	mu->resize(num_ang_store);
	eta->resize(num_ang_store);
	d_wt.resize(num_ang_store);

	// All angles have a constant weight
	double wt = four_pi / num_ang_store;

	// Octant directions
	double s[8][3] = { {1,1,1},{-1,1,1},{1,-1,1},{-1,-1,1},
			           {1,1,-1},{-1,1,-1},{1,-1,-1},{-1,-1,-1} };

	// Fill the quadrature angle and weight arrays
	int k = 0;
	for ( int o = 0; o < num_oct; o++ )
	{
	    for ( int i = 0; i < num_xi_2; i++ )
	    {
	        for ( int j = 0; j < num_omega_4; j++ )
	        {
	            (*mu)[k]  	= s[o][MU]*tmp_mu[j][i];
	            (*eta)[k] 	= s[o][ETA]*tmp_eta[j][i];
	            d_wt[k] 	= wt;
	            ++k;
	        }
	    }
	}
	Ensure ( k == num_ang_store );

	// Normalize the quadrature set
	double wsum = 0.0;
	for (size_t n = 0; n <= num_ang_store - 1; ++n)
		wsum = wsum + d_wt[n];

	for (size_t n = 0; n <= num_ang_store - 1; ++n)
		d_wt[n] = d_wt[n] * (norm / wsum);

#ifdef DEBUG
	display();
#endif

	Ensure (d_angles[MU].size() == num_ang_store);
	Ensure (d_angles[ETA].size() == num_ang_store);
}

//---------------------------------------------------------------------------//
// PUBLIC FUNCTIONS
//---------------------------------------------------------------------------//
/*!
 * \brief Display the quadrature.
 */
//void Uniform_Equal_3D::display() const
//{
//
//	using std::cout;
//	using std::endl;
//	using std::printf;
//	Require ( (d_angles.size()==2) || (d_angles.size()==3) );
//	cout << endl << "The Quadrature directions and weights are:"
//			<< endl << endl;
//	if ( d_angles.size() == 2 )
//	{
//		cout << "   m         mu              eta             wt       " << endl;
//		cout << "  ---   --------------  --------------  --------------" << endl;
//	}
//	else if ( d_angles.size() == 3 )
//	{
//		cout << "   m         mu              eta             xi              wt       " << endl;
//		cout << "  ---   --------------  --------------  --------------  --------------" << endl;
//	}
//	double sum_wt = 0.0;
//	for ( int ix = 0; ix < d_angles[MU].size(); ++ix )
//	{
//		if ( d_angles.size() == 2 )
//		{
//			printf ("%4i    %13.10f   %13.10f   %13.10f   \n",
//					ix, d_angles[MU][ix], d_angles[ETA][ix], d_wt[ix] );
//		}
//		else
//		{
//			printf ("%4i    %13.10f   %13.10f   %13.10f   %13.10f   \n",
//					ix, d_angles[MU][ix], d_angles[ETA][ix], d_angles[XI][ix], d_wt[ix] );
//		}
//		sum_wt += d_wt[ix];
//	}
//	cout << endl << "  The sum of the weights is " << sum_wt << endl;
//	cout << endl;
//
//}

} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of Uniform_Equal_3D.cc
//---------------------------------------------------------------------------//
